1 Setup

Libraries

suppressPackageStartupMessages({
  library(DESeq2)
  library(dplyr)
  library(gplots)
  library(ggplot2)
  library(here)
  library(hyperSpec)
  library(limma)
  library(readr)
  library(tibble)
  library(RColorBrewer)
  library(VennDiagram)
})

Helpers

source(here("UPSCb-common/src/R/gopher.R"))

Graphics

hpal <- colorRampPalette(c("blue","white","red"))(100)
pal <- brewer.pal(8,"Dark2")

Data

suppressPackageStartupMessages(load(here("analysis/batchCorrection/vsd.rda")))

Time points of interest

sample.sel <- (vsd$Time %in% paste0("B",3:6) & vsd$Tissue =="SE") | 
(vsd$Tissue!="SE" & vsd$Time %in% paste0("B",4:9))
vsd <- vsd[,sample.sel]
vsd$Time <- droplevels(vsd$Time)

Devise a pseudotime from the PCA after bacth correction

# The order is B4-5 ZE; B3 SE; B6 ZE, B4-5 SE, B7-9 ZE, B6 SE
# corresponding to pseudo-time 1-2,3,4,5-6,7-9,10
vsd$PseudoTime <- as.integer(vsd$Time) -1 + 
  ifelse(vsd$Tissue=="SE",
         ifelse(vsd$Time=="B3",3,
                ifelse(vsd$Time %in% c("B4","B5"),4,7)),
         ifelse(vsd$Time=="B6",1,
           ifelse(vsd$Time %in% paste0("B",7:9),3,0)))

Gene of interests

goi <- read_tsv(here("doc/MA_list_known_SE_genes.txt"),
                col_types=cols(.default=col_character()))

2 Functions

"line_plot" <- function(vsd=vsd,gene_id=gene_id,gene_name=gene_name){
  message(paste("Plotting",gene_id))
  sel <- grepl(gene_id,rownames(vsd))
  stopifnot(sum(sel)==1)
  
  p <- ggplot(bind_cols(as.data.frame(colData(vsd)),
                        data.frame(value=assay(vsd)[sel,])),
              aes(x=PseudoTime,y=value,col=Tissue,group=Tissue)) +
    geom_point() + geom_smooth() +
    scale_y_continuous(name="VST expression") + 
    ggtitle(label=paste("Expression for: ",gene_id,"(",gene_name,")"))
  
  suppressMessages(suppressWarnings(plot(p)))
  return(NULL)
}

3 Plots

dev.null <- apply(goi,1,function(ro){line_plot(vsd,ro[2],ro[1])})

## Plotting MA_86195g0010

## Plotting MA_98095g0010

## Plotting MA_121578g0010

## Plotting MA_903853g0010

## Plotting MA_52791g0010

## Plotting MA_216909g0010

## Plotting MA_68722g0010

## Plotting MA_147422g0010

## Plotting MA_823242g0010

## Plotting MA_66505g0010

## Plotting MA_66505g0020

## Plotting MA_4929831g0010

## Plotting MA_104203g0010

## Plotting MA_10428962g0010

## Plotting MA_196219g0010

## Plotting MA_8832398g0010

## Plotting MA_42637g0010

## Plotting MA_8990972g0010

## Plotting MA_8227330g0010

## Plotting MA_10429361g0010

## Plotting MA_10434655g0010

## Plotting MA_10432909g0020

## Plotting MA_119006g0010

## Plotting MA_95218g0010

## Plotting MA_10431854g0020

## Plotting MA_255g0010

## Plotting MA_10429533g0010

## Plotting MA_113269g0010

## Plotting MA_166600g0010

## Plotting MA_10163020g0010

## Plotting MA_69685g0020

## Plotting MA_197296g0010

## Plotting MA_58651g0010

## Plotting MA_10435775g0010

## Plotting MA_141759g0010

## Plotting MA_15508g0010

## Plotting MA_9199701g0010

## Plotting MA_39847g0010

## Plotting MA_10428957g0010

## Plotting MA_10436994g0010

## Plotting MA_117784g0010

## Plotting MA_808776g0010

## Plotting MA_8273259g0010

## Plotting MA_10087005g0010

## Plotting MA_10431301g0020

## Plotting MA_176924g0010

## Plotting MA_623831g0010

## Plotting MA_12471g0030

## Plotting MA_102140g0010

## Plotting MA_9725777g0030

## Plotting MA_10429954g0020

## Plotting MA_10437050g0010

## Plotting MA_135287g0010

## Plotting MA_10430507g0030

## Plotting MA_10360330g0010

## Plotting MA_10386237g0010

## Plotting MA_9968141g0010

## Plotting MA_35079g0010

## Plotting MA_1036215g0010

## Plotting MA_10053571g0010

## Plotting MA_1090407g0010

## Plotting MA_17729g0010

## Plotting MA_10434917g0010

## Plotting MA_10437070g0010

## Plotting MA_10435722g0010

## Plotting MA_10435722g0010

## Plotting MA_66505g0010

## Plotting MA_66505g0020

4 Differential expression

We assume the model to be expression as a function of Stage and Tissue and their interaction i.e. the two variables are not considered independent. However, we will group the Stages, to compare Maturation (SE B3-5 and ZE B3-6) and Desiccation (SE B6 and ZE B7-9).

vsd$Stage<-ifelse(vsd$PseudoTime<7,"Maturation","Dessication")
Stage <- vsd$Stage
Tissue<-vsd$Tissue
design <- model.matrix(~Stage*Tissue)

remove genes with t0o little variance

mat <- assay(vsd)[rowMads(assay(vsd))>0,]

fit <- lmFit(mat, design)
fit <- eBayes(fit)

The possible coefficients

colnames(design)

## [1] "(Intercept)"              "StageMaturation"         
## [3] "TissueSE"                 "TissueZE"                
## [5] "StageMaturation:TissueSE" "StageMaturation:TissueZE"

4.1 ZE vs. SE (maturation)

fit2 <- contrasts.fit(fit, c(0,0,0,0,-1,1))
fit2 <- eBayes(fit2)
ZSM <- topTable(fit2,p.value=0.01,lfc=0.5,adjust.method="BH",n=nrow(vsd))
sample.sel <- vsd$Stage == "Maturation" & vsd$Tissue != "FMG"
heatmap.2(t(scale(t(mat[rownames(ZSM),sample.sel]))),
          distfun=pearson.dist,
          hclustfun=function(X){hclust(X,method="ward.D2")},
          labRow = NA,trace = "none",
          labCol = paste(vsd$Tissue,vsd$Time)[sample.sel],
          col=hpal,cexCol=0.8)

4.2 ZE vs. SE (desiccation)

fit2 <- contrasts.fit(fit, c(0,0,-1,1,0,0))
fit2 <- eBayes(fit2)
ZSD <- topTable(fit2,p.value=0.01,lfc=0.5,adjust.method="BH",n=nrow(vsd))
sample.sel <- vsd$Stage != "Maturation" & vsd$Tissue != "FMG"
heatmap.2(t(scale(t(mat[rownames(ZSD),sample.sel]))),
          distfun=pearson.dist,
          hclustfun=function(X){hclust(X,method="ward.D2")},
          labRow = NA,trace = "none",
          labCol = paste(vsd$Tissue,vsd$Time)[sample.sel],
          col=hpal,cexCol=0.8)

4.3 Overlap

grid.newpage()
grid.draw(venn.diagram(list(ZSD=rownames(ZSD),
                            ZSM=rownames(ZSM)),NULL,fill=pal[1:2]))

5 Gene Ontology Enrichment

GO <- lapply(lapply(list(rownames(ZSD),rownames(ZSM)),
                    sub,pattern="\\.1$",replacement=""),
             gopher,background=sub("\\.1","",rownames(mat)),
             task=list("go"),alpha=0.01,
             url="pabies")

## Loading required package: tidyr

## 
## Attaching package: 'tidyr'

## The following object is masked from 'package:S4Vectors':
## 
##     expand

## Loading required package: jsonlite

6 Export

dir.create(here("analysis/limma"),showWarnings=FALSE)  
write_delim(ZSD %>% rownames_to_column("ID"),here("analysis/limma/ZE-vs-SE-desiccation_DE-results.tsv"))
write_delim(ZSM %>% rownames_to_column("ID"),here("analysis/limma/ZE-vs-SE-maturation_DE-results.tsv"))
write_delim(GO[[1]]$go,here("analysis/limma/ZE-vs-SE-desiccation_DE_GO-results.tsv"))
write_delim(GO[[2]]$go,here("analysis/limma/ZE-vs-SE-maturation_DE_GO-results.tsv"))
write_delim(GO[[1]]$go[,c("id","padj")],here("analysis/limma/ZE-vs-SE-desiccation_DE_GO-for-REVIGO-results.tsv"))
write_delim(GO[[2]]$go[,c("id","padj")],here("analysis/limma/ZE-vs-SE-maturation_DE_GO-for-REVIGO-results.tsv"))

7 SessionInfo

## R version 4.0.3 (2020-10-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.5 LTS
## 
## Matrix products: default
## BLAS/LAPACK: /usr/lib/x86_64-linux-gnu/libopenblasp-r0.2.20.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
##  [1] grid      parallel  stats4    stats     graphics  grDevices utils    
##  [8] datasets  methods   base     
## 
## other attached packages:
##  [1] jsonlite_1.7.2              tidyr_1.1.2                
##  [3] VennDiagram_1.6.20          futile.logger_1.4.3        
##  [5] RColorBrewer_1.1-2          tibble_3.0.4               
##  [7] readr_1.4.0                 limma_3.46.0               
##  [9] hyperSpec_0.99-20201127     xml2_1.3.2                 
## [11] lattice_0.20-41             here_1.0.1                 
## [13] ggplot2_3.3.3               gplots_3.1.1               
## [15] dplyr_1.0.2                 DESeq2_1.30.0              
## [17] SummarizedExperiment_1.20.0 Biobase_2.50.0             
## [19] MatrixGenerics_1.2.0        matrixStats_0.57.0         
## [21] GenomicRanges_1.42.0        GenomeInfoDb_1.26.2        
## [23] IRanges_2.24.1              S4Vectors_0.28.1           
## [25] BiocGenerics_0.36.0        
## 
## loaded via a namespace (and not attached):
##  [1] nlme_3.1-151           bitops_1.0-6           bit64_4.0.5           
##  [4] httr_1.4.2             rprojroot_2.0.2        tools_4.0.3           
##  [7] R6_2.5.0               KernSmooth_2.23-18     DBI_1.1.0             
## [10] lazyeval_0.2.2         mgcv_1.8-33            colorspace_2.0-0      
## [13] withr_2.3.0            tidyselect_1.1.0       curl_4.3              
## [16] bit_4.0.4              compiler_4.0.3         formatR_1.7           
## [19] DelayedArray_0.16.0    labeling_0.4.2         caTools_1.18.0        
## [22] scales_1.1.1           genefilter_1.72.0      stringr_1.4.0         
## [25] digest_0.6.27          rmarkdown_2.6          XVector_0.30.0        
## [28] jpeg_0.1-8.1           pkgconfig_2.0.3        htmltools_0.5.0       
## [31] highr_0.8              rlang_0.4.10           RSQLite_2.2.1         
## [34] farver_2.0.3           generics_0.1.0         BiocParallel_1.24.1   
## [37] gtools_3.8.2           RCurl_1.98-1.2         magrittr_2.0.1        
## [40] GenomeInfoDbData_1.2.4 Matrix_1.3-0           Rcpp_1.0.5            
## [43] munsell_0.5.0          lifecycle_0.2.0        stringi_1.5.3         
## [46] yaml_2.2.1             zlibbioc_1.36.0        blob_1.2.1            
## [49] crayon_1.3.4           splines_4.0.3          annotate_1.68.0       
## [52] hms_0.5.3              locfit_1.5-9.4         knitr_1.30            
## [55] pillar_1.4.7           geneplotter_1.68.0     futile.options_1.0.1  
## [58] XML_3.99-0.5           glue_1.4.2             evaluate_0.14         
## [61] latticeExtra_0.6-29    lambda.r_1.2.4         vctrs_0.3.6           
## [64] png_0.1-7              testthat_3.0.1         gtable_0.3.0          
## [67] purrr_0.3.4            xfun_0.20              xtable_1.8-4          
## [70] survival_3.2-7         AnnotationDbi_1.52.0   memoise_1.1.0         
## [73] ellipsis_0.3.1

ZE and SE differential expression analysis

Nicolas Delhomme && Michael Stewart

2021-01-14